{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbe99385",
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m内核无法启动，因为 Python 环境“env4 (Python -1.-1.-1)”不再可用。请考虑选择另一个内核或刷新 Python 环境列表。"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Point\n",
    "from shapely.ops import unary_union\n",
    "import pickle\n",
    "import time\n",
    "import traceback\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "try:\n",
    "    from pso_core import HierarchicalRoadNetwork, OptimizedMultiSitePSO, create_scenario_map\n",
    "    print(\"✅ 核心模块 'pso_core.py' 导入成功。\")\n",
    "except ImportError:\n",
    "    print(\"❌ 错误：无法从 'pso_core.py' 导入核心类。\")\n",
    "    print(\"   请确保您的所有类（HierarchicalRoadNetwork, OptimizedMultiSitePSO, etc.）都在一个名为 'pso_core.py' 的文件中，\")\n",
    "    print(\"   并且这个文件与当前的Notebook在同一个目录下。\")\n",
    "\n",
    "# 配置Matplotlib以正确显示中文和负号\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']\n",
    "plt.rcParams['axes.unicode_minus'] = False\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "976ac9a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    print(\"--- 开始加载和预处理所有地理数据 ---\")\n",
    "    \n",
    "    # 1. 导入坐标转换库\n",
    "    from coord_convert.transform import gcj2wgs\n",
    "\n",
    "    # 2. 定义文件路径\n",
    "    boundary_shapefile = \"仓山区地图边界.shp\"\n",
    "    poi_csv = \"poi仓山区.csv\"\n",
    "    population_shapefile = \"population.shp\"\n",
    "    ruixing_csv = \"福州市仓山区瑞幸咖啡坐标.csv\"\n",
    "    road_shapefile = \"2025年福建省道路数据.shp\"\n",
    "    road_cache_path = \"road_cache.pkl\"\n",
    "\n",
    "    # 3. 加载并处理【边界】数据\n",
    "    boundary_gdf_raw = gpd.read_file(boundary_shapefile)\n",
    "    if boundary_gdf_raw.crs is None:\n",
    "        boundary_gdf_raw.set_crs(\"EPSG:4326\", inplace=True)\n",
    "    boundary_gdf_3857 = boundary_gdf_raw.to_crs(\"EPSG:3857\")\n",
    "    boundary_polygon_3857 = boundary_gdf_3857.union_all()\n",
    "    print(\"✅ 边界数据处理完毕。\")\n",
    "\n",
    "    # 4. 加载并处理【POI】数据\n",
    "    print(\"正在转换POI坐标...\")\n",
    "    poi_df = pd.read_csv(poi_csv, encoding='gbk')\n",
    "    corrected_poi_coords = [gcj2wgs(lon, lat) for lon, lat in zip(poi_df.geolocation_lng, poi_df.geolocation_lat)]\n",
    "    poi_gdf_3857 = gpd.GeoDataFrame(\n",
    "        geometry=[Point(lon, lat) for lon, lat in corrected_poi_coords],\n",
    "        crs=\"EPSG:4326\"\n",
    "    ).to_crs(\"EPSG:3857\")\n",
    "    poi_coords_3857 = [(p.x, p.y) for p in poi_gdf_3857.geometry]\n",
    "    print(\"✅ POI数据处理完毕。\")\n",
    "\n",
    "    # 5. 加载并处理【人口】数据\n",
    "    population_gdf = gpd.read_file(population_shapefile)\n",
    "    population_gdf_3857 = population_gdf.to_crs(\"EPSG:3857\")\n",
    "    population_points_3857 = [\n",
    "        ((point.x, point.y), value) \n",
    "        for point, value in zip(population_gdf_3857.geometry, population_gdf_3857['population'])\n",
    "        if value > 0\n",
    "    ]\n",
    "    print(f\"✅ 人口数据处理完毕，共加载 {len(population_points_3857)} 个有效人口点。\")\n",
    "\n",
    "    # 6. 加载并【修正】瑞幸咖啡坐标\n",
    "    print(\"正在转换瑞幸咖啡坐标 (GCJ-02 -> WGS-84)...\")\n",
    "    ruixing_df = pd.read_csv(ruixing_csv, encoding='gbk')\n",
    "    corrected_ruixing_coords = [gcj2wgs(lon, lat) for lon, lat in zip(ruixing_df.geolocation_lng, ruixing_df.geolocation_lat)]\n",
    "    ruixing_corrected_lng = [coord[0] for coord in corrected_ruixing_coords]\n",
    "    ruixing_corrected_lat = [coord[1] for coord in corrected_ruixing_coords]\n",
    "    ruixing_gdf_wgs84 = gpd.GeoDataFrame(\n",
    "        ruixing_df, \n",
    "        geometry=gpd.points_from_xy(ruixing_corrected_lng, ruixing_corrected_lat),\n",
    "        crs=\"EPSG:4326\"\n",
    "    )\n",
    "    print(\"✅ 瑞幸咖啡坐标加载并修正成功。\")\n",
    "\n",
    "    # --- 7. 加载和处理【路网】数据 (【关键修正】) ---\n",
    "    print(\"正在处理路网数据...\")\n",
    "    boundary_gdf_wgs84 = boundary_gdf_raw.to_crs(\"EPSG:4326\") # 确保用于裁剪的边界是WGS84\n",
    "    road_gdf_raw = gpd.read_file(road_shapefile, bbox=boundary_gdf_wgs84)\n",
    "    road_gdf_clipped = gpd.clip(road_gdf_raw, boundary_gdf_wgs84)\n",
    "    \n",
    "    # a. 先将整个 GeoDataFrame 转换为目标坐标系\n",
    "    road_gdf_clipped_3857 = road_gdf_clipped.to_crs(\"EPSG:3857\")\n",
    "    \n",
    "    # b. 创建一个副本以避免 SettingWithCopyWarning\n",
    "    road_gdf_simplified_3857 = road_gdf_clipped_3857.copy()\n",
    "    \n",
    "    # c. 🔥 在副本的 'geometry' 列上进行 simplify 操作\n",
    "    # 这种方式可以确保 road_gdf_simplified_3857 变量本身始终是一个 GeoDataFrame\n",
    "    road_gdf_simplified_3857['geometry'] = road_gdf_clipped_3857.geometry.simplify(tolerance=10)\n",
    "    \n",
    "    print(\"✅ 路网数据处理完毕。\")\n",
    "    \n",
    "    # 8. 实例化路网对象\n",
    "    print(\"\\n--- 正在实例化路网对象并生成/更新缓存 ---\")\n",
    "    road_network = HierarchicalRoadNetwork(road_gdf_simplified_3857, cache_file=road_cache_path)\n",
    "    print(\"✅ 路网对象实例化完成。\")\n",
    "\n",
    "    # 9. 打包并保存输入数据\n",
    "    pso_input_data = {\n",
    "        \"boundary_polygon_3857\": boundary_polygon_3857,\n",
    "        \"poi_coords_3857\": poi_coords_3857,\n",
    "        \"population_points_3857\": population_points_3857,\n",
    "        \"road_network_cache\": road_cache_path \n",
    "    }\n",
    "    with open(\"pso_input_data.pkl\", \"wb\") as f:\n",
    "        pickle.dump(pso_input_data, f)\n",
    "\n",
    "    print(\"\\n\\n🎉 [成功] 所有输入数据已成功打包并保存到 'pso_input_data.pkl'\")\n",
    "\n",
    "except Exception as e:\n",
    "    print(f\"\\n❌ [严重错误] 在数据准备阶段发生错误: {e}\")\n",
    "    traceback.print_exc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30463b7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import time\n",
    "\n",
    "print(\"=\"*80)\n",
    "print(\"🚀 即将开始在外部Python进程中执行多进程PSO优化...\")\n",
    "print(\"💻 脚本的所有输出（包括日志和错误）都将实时显示在下方。\")\n",
    "print(\"=\"*80)\n",
    "\n",
    "# --- 1. 定义要执行的命令 ---\n",
    "command = [\"python\", \"-u\", \"run_pso_multiprocess.py\"]\n",
    "# \"-u\" 参数是关键：它强制Python使用无缓冲的输出流，\n",
    "# 确保外部脚本的每一次 print() 或 logger 输出都立刻被我们捕获，而不是等待缓冲区满。\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "# --- 2. 使用 subprocess.Popen 实时捕获和打印输出 ---\n",
    "try:\n",
    "    # 启动子进程，并将stdout和stderr合并\n",
    "    process = subprocess.Popen(\n",
    "        command,\n",
    "        stdout=subprocess.PIPE,\n",
    "        stderr=subprocess.STDOUT, # 🔥 关键修改：将错误流合并到标准输出流\n",
    "        text=True,\n",
    "        encoding='utf-8',\n",
    "        errors='replace' # 防止因编码问题中断\n",
    "    )\n",
    "\n",
    "    # 实时读取并打印子进程的输出\n",
    "    while True:\n",
    "        output = process.stdout.readline()\n",
    "        if output == '' and process.poll() is not None:\n",
    "            break\n",
    "        if output:\n",
    "            print(output.strip()) # 实时打印每一行输出\n",
    "\n",
    "    # 等待进程结束并获取返回码\n",
    "    return_code = process.poll()\n",
    "\n",
    "finally:\n",
    "    end_time = time.time()\n",
    "    total_minutes = (end_time - start_time) / 60\n",
    "    print(\"\\n\" + \"=\"*80)\n",
    "    if return_code == 0:\n",
    "        print(f\"🎉 [成功] 外部脚本执行完毕。总耗时: {total_minutes:.2f} 分钟。\")\n",
    "        print(\"   优化结果已保存到 'pso_results.pkl' 文件中。\")\n",
    "    else:\n",
    "        print(f\"❌ [错误] 外部脚本执行失败！返回码: {return_code}\")\n",
    "        print(\"   请仔细阅读上面的输出日志，找到导致脚本崩溃的错误信息。\")\n",
    "    print(\"=\"*80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d14599cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# =========================================================================\n",
    "# === 单元格 4: 加载、分析与可视化结果 (【最终整合版】) ====================\n",
    "# =========================================================================\n",
    "# 描述:\n",
    "# 这个单元格负责加载由后台脚本生成的优化结果，进行深入的量化分析，\n",
    "# 并生成两个核心的可视化产出：\n",
    "# 1. PSO算法的两阶段收敛曲线图。\n",
    "# 2. 包含最优选址点、现有门店、人口热力等的最终可交互地图。\n",
    "# =========================================================================\n",
    "\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely.geometry import Point\n",
    "import traceback\n",
    "import os\n",
    "\n",
    "# --- 假设以下变量和函数已在前面的单元格中定义并可访问 ---\n",
    "# ruixing_csv (瑞幸门店CSV文件路径)\n",
    "# gcj2wgs (坐标转换函数)\n",
    "# create_scenario_map (地图生成函数)\n",
    "# boundary_polygon_3857 (区域边界)\n",
    "# poi_coords_3857 (POI坐标)\n",
    "# population_points_3857 (人口点)\n",
    "# road_gdf_simplified_3857 (路网GeoDataFrame)\n",
    "\n",
    "print(\"\\n--- 开始加载和分析优化结果 ---\")\n",
    "\n",
    "try:\n",
    "    with open(\"pso_results.pkl\", \"rb\") as f:\n",
    "        results = pickle.load(f)\n",
    "    print(\"✅ 结果文件 'pso_results.pkl' 加载成功。\")\n",
    "\n",
    "    # 2. 从结果字典中解包数据\n",
    "    optimal_locations_3857 = results[\"optimal_locations_3857\"]\n",
    "    analysis_metrics = results[\"analysis_metrics\"]\n",
    "    convergence_history = results[\"convergence_history\"] # 现在这是一个简单的数值列表\n",
    "\n",
    "    # 3. 将米制坐标转换为WGS84经纬度\n",
    "    optimal_gdf_3857 = gpd.GeoDataFrame(geometry=[Point(loc) for loc in optimal_locations_3857], crs=\"EPSG:3857\")\n",
    "    optimal_gdf_wgs84 = optimal_gdf_3857.to_crs(\"EPSG:4326\")\n",
    "    optimal_locations_wgs84 = [(p.x, p.y) for p in optimal_gdf_wgs84.geometry]\n",
    "\n",
    "    # 4. 打印一份清晰、详细的最终报告\n",
    "    print(\"\\n\" + \"=\"*40 + \" ✨ 最终方案分析报告 ✨ \" + \"=\"*40)\n",
    "    print(f\"\\n--- 优化摘要 ---\")\n",
    "    print(f\"  - 最终找到 {len(optimal_locations_wgs84)} 个最优站点位置。\")\n",
    "    model_score = analysis_metrics.get('model_score', {})\n",
    "    print(f\"  - 模型最终得分 (越低越好): {model_score.get('total', 'N/A'):.4f}\")\n",
    "    \n",
    "    print(\"\\n--- 关键业务指标 ---\")\n",
    "    pop_coverage = analysis_metrics.get('population_coverage', {})\n",
    "    poi_service = analysis_metrics.get('poi_service', {})\n",
    "    spacing_stats = analysis_metrics.get('spacing_stats', {})\n",
    "    print(f\"  - 📈 实际人口覆盖率: {pop_coverage.get('coverage_ratio', 0) * 100:.2f}%\")\n",
    "    print(f\"  - 🎯 POI服务率: {poi_service.get('service_ratio', 0) * 100:.2f}%\")\n",
    "    print(f\"  - 📏 站点平均间距: {spacing_stats.get('average_distance_m', 0):.0f} 米\")\n",
    "    print(f\"  - ⚠️ 站点最小间距: {spacing_stats.get('min_distance_m', 0):.0f} 米\")\n",
    "\n",
    "    print(\"\\n--- 成本构成 (0-1, 越低越好) ---\")\n",
    "    print(f\"  - POI便利性成本 (加权):   {model_score.get('poi_weighted', 'N/A'):.4f}\")\n",
    "    print(f\"  - 人口覆盖成本 (加权):     {model_score.get('population_weighted', 'N/A'):.4f}\")\n",
    "    print(f\"  - 站点间距成本 (加权):     {model_score.get('spacing_weighted', 'N/A'):.4f}\")\n",
    "    \n",
    "    print(\"\\n--- 最优站点坐标 (WGS84 经纬度) ---\")\n",
    "    for i, (lon, lat) in enumerate(optimal_locations_wgs84):\n",
    "        print(f\"  - 站点 {i+1:02d}: (经度 {lon:.6f}, 纬度 {lat:.6f})\")\n",
    "    print(\"=\"*100)\n",
    "    \n",
    "    # 5. 生成并显示两个阶段的收敛曲线图\n",
    "    print(\"\\n正在生成收敛曲线图...\")\n",
    "    if not convergence_history or not all(isinstance(x, (int, float)) for x in convergence_history):\n",
    "        print(\"⚠️ 未能生成收敛曲线图，因为历史分数数据格式不正确或为空。\")\n",
    "    else:\n",
    "        plt.figure(figsize=(10, 6))\n",
    "        \n",
    "        # 🔥 直接绘制收敛历史，不再需要处理阶段切换\n",
    "        plt.plot(convergence_history, marker='o', markersize=3, linewidth=1.5)\n",
    "        \n",
    "        plt.title(\"PSO 优化收敛曲线 (高精度模式)\", fontsize=16) # 标题可以更新\n",
    "        plt.xlabel(\"迭代次数\", fontsize=12)\n",
    "        plt.ylabel(\"目标函数值 (成本)\", fontsize=12)\n",
    "        plt.grid(True, linestyle='--', alpha=0.6)\n",
    "        plt.show()\n",
    "\n",
    "    \n",
    "    print(\"\\n正在生成最终结果的可交互地图...\")\n",
    "    ruixing_df = pd.read_csv(ruixing_csv, encoding='gbk')\n",
    "    corrected_ruixing_coords = [gcj2wgs(lon, lat) for lon, lat in zip(ruixing_df.geolocation_lng, ruixing_df.geolocation_lat)]\n",
    "    corrected_lng = [coord[0] for coord in corrected_ruixing_coords]\n",
    "    corrected_lat = [coord[1] for coord in corrected_ruixing_coords]\n",
    "    ruixing_gdf_wgs84 = gpd.GeoDataFrame(\n",
    "        ruixing_df, \n",
    "        geometry=gpd.points_from_xy(corrected_lng, corrected_lat),\n",
    "        crs=\"EPSG:4326\"\n",
    "    )\n",
    "    \n",
    "    create_scenario_map(\n",
    "        scenario_name=\"PSO_MultiProcess_Final\",\n",
    "        optimal_locations=optimal_locations_3857,\n",
    "        boundary_polygon_3857=boundary_polygon_3857,\n",
    "        poi_coords_3857=poi_coords_3857,\n",
    "        population_points_3857=population_points_3857,\n",
    "        road_network_gdf_3857=road_gdf_simplified_3857,\n",
    "        actual_locations_gdf=ruixing_gdf_wgs84,\n",
    "        output_folder=\"final_maps\"\n",
    "    )\n",
    "    print(\"\\n🎉 [完成] 所有分析和可视化任务已结束。\")\n",
    "\n",
    "except FileNotFoundError:\n",
    "    print(\"\\n❌ 错误：找不到结果文件 'pso_results.pkl'。\")\n",
    "    print(\"   请确认您已经成功运行了【单元格 3】并且没有发生错误。\")\n",
    "except Exception as e:\n",
    "    print(f\"\\n❌ [严重错误] 在加载或分析结果时发生错误: {e}\")\n",
    "    traceback.print_exc()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env4",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "-1.-1.-1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
